home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Precision Software Appli…tions Silver Collection 4
/
Precision Software Applications Silver Collection Volume 4 (1993).iso
/
stats
/
chadyn.exe
/
YARON.C
< prev
next >
Wrap
C/C++ Source or Header
|
1988-11-28
|
2KB
|
90 lines
/********************** YDblVdp for double vanderPol *************************/
/********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
#include "yinclud.h"
aron(k, Y) /* 2 coupled vanderpol oscillators -- see also map 23 */
double k[],
Y[];
{
double bb,
tt,
t,
X,
W,
yy,
N,
exp(), expY, expXYW, expX, expWY, expN;
t = Y[0];
X = Y[1];
W = Y[2];
yy = Y[3];
N = Y[4];
bb = C2;
tt = 1.+ beta * cos(C4 * t);
expY = -(rho - 1.) * C1 * tt * exp(yy);
expXYW = (C3 + C1) * tt * exp(X + yy - W);
expX = C1 * rho * bb * exp(-X);
expWY = (C5 + C1) * exp(W - yy);
expN = C1 * bb * exp(-N);
k[0] = 1.;
k[1] = expY + expX - C1;
k[2] = expXYW - (C3 + C1);
k[3] = expWY - (C5 + C1);
k[4] = expN - C1;
#ifdef GARBAGE
if(num_lyap > 0) {
for(i = 0; i < num_lyap; i++) {
base = lyapzero + vec_dim * i;
k[base] = Y[base + 1];
k[base + 1] = () * Y[base]
+ C1 * (1.- X1 * X1) * Y[base + 1]
+ rho * Y[base + 2];
k[base + 2] = Y[base + 3];
k[base + 3] =
rho * Y[base]
+ (-2 * C5 * X2 * Y2 - C6 - 3 * C7 * X2 * X2) * Y[base + 2]
+ C5 * (1 - X2 * X2) * Y[base + 3];
}
}
#endif
}
init_aron() {
extern int aron();
zeroth = 1;
vec_dim = 3; /* the dimension of the Lyapunov vectors =
phase space dim; the component y[4] is not
being counted in this because it drifts
slowly */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 5; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
dim = lyapzero + num_lyap * vec_dim;
/* needed for rungekutta */
X_upper =.3; /* x scale */
X_lower = -.3;
Y_upper = 1.5; /* y scale */
Y_lower = -5.;
X_coord = 1;
Y_coord = 3;
C1 =.02;
C2 =.6667;
C3 = 35.84;
C4 = twopi;
C5 = 100.;
beta =.1;
rho = 18.;
steps_per_cycle = 200;
step = (twopi / C4) / steps_per_cycle;
DE = aron; /* DE is a pointer to a function */
y[eqn0 + 1] = -0.08960;
y[eqn0 + 2] = -1.90847;
y[eqn0 + 3] = -1.90918;
y[eqn0 + 4] = -0.00662;
setequal(0, 1, eqn1); /* reinitializes y[] */
}